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Abstract In vitro cell biology assays play a crucial role in informing our un- 
derstanding of the migratory, proliferative and invasive properties of many 
cell types in different biological contexts. While mono-culture assays involve 
the study of a population of cells composed of a single cell type, co-culture 
assays study a population of cells composed of multiple cell types (or subpop- 
ulations of cells). Such co-culture assays can provide more realistic insights 
into many biological processes including tissue repair, tissue regeneration and 
malignant spreading. Typically, system parameters, such as motility and pro- 
liferation rates, are estimated by calibrating a mathematical or computational 
model to the observed experimental data. However, parameter estimates can 
be incredibly sensitive to the choice of model and modelling framework. This 
observation motivates us to consider the fundamental question of how we can 
best choose a model to facilitate accurate parameter estimation for a particu- 
lar assay. In this work we describe three mathematical models of mono-culture 
and co-culture assays that include different levels of spatial detail. We study 
various spatial summary statistics to explore if they can be used to distin- 
guish between the suitability of each model over a range of parameter space. 
Our results for mono-culture experiments are promising, in that we suggest 
two spatial statistics that can be used to direct model choice. However, co- 
culturc experiments are far more challenging: we show that these same spatial 
statistics which provide useful insight into mono-culture systems are insuffi- 
cient for co-culture systems. Therefore, we conclude that great care ought to 
be exercised when estimating the parameters of co-culture assays. 
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Fig. 1 Snapshots showing two growth-to-confluencc assays. Images in (a)— (c) show a mono- 
culture experiment using a population of 3T3 fibroblast cells whereas the images in (d)-(f) 
show a co-culture experiment containing a subpopulation of melanocyte cells (red arrows) 
amongst a population of keratinocytes (no arrows). The scale bar in all subfigures corre- 
sponds to 100 fim. 



1 Introduction 

In vitro cell biology assays are an essential element in the study of the migra- 
tory, proliferative and invasive properties of different types of cells [8] , and they 
provide insight into various phenomena including malignant spreading [25] and 
wound healing [27]. While many types of such in vitro assays involve a mono- 
culture system of a single cell type, many other applications require the anal- 
ysis of a co-culture system that investigates the migration and proliferation of 
multiple cell types or subpopulations. For example, wound healing requires the 
controlled proliferation and migration of both keratinocytes and fibroblasts as 
well as a number of complex interactions between these two cell types [26] . 

An example of an in vitro mono-culture assay is shown in Figure l(a)-(c). 
This is a growth-to-confluence assay involving 3T3 fibroblast cells [22]. In this 
kind of assay, an initially uniform population of cells is placed on a tissue 
culture plate and monitored in real time as the individual cells within the 
population move and proliferate to eventually form a confluent monolayer [17]. 
Analyzing the rate at which the density of cells increases with time allows us 
to make inferences about the rate at which the cells proliferate [17], which is 
an essential component of collective cell spreading. Alternatively, an example 
of an in vitro co-culture assay is shown in Figure l(d)-(f). This is a growth-to- 
confluence assay involving two cell types: melanocytes and keratinocytes. In 
this assay, the two cell types are placed on a tissue culture plate and monitored 
in real time as cells from both subpopulations move and proliferate, eventually 
leading to a confluent monolayer of cells containing a mixture of both cell types. 
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Our previous work has focused on developing and applying mathematical 
models to interpret mono-culture growth-to-confluence experiments [10,17]. 
In particular, we showed that mono-culture growth-to-confluence experiments 
can be described using three different mathematical modelling frameworks. 
Firstly, we considered a stochastic description of individual cell motility, pro- 
liferation and death events which has the advantage of directly incorporating 
individual cell-level behaviours and naturally gives rise to spatial correlations 
in cell locations, but is analytically intractable [4,5]. Secondly we considered 
the traditional corresponding mean-field description of the average cell pop- 
ulation density which has the advantage of being analytically tractable but 
suffers from the disadvantage of neglecting spatial correlations in the distri- 
bution of individual cells [7,9]. Thirdly, we considered the more sophisticated 
corresponding moment dynamics description of the average cell density which 
has the advantage of being computationally tractable and approximately in- 
corporating the effects of spatial correlations in the distribution of individual 
cells in the population [11]. 

Previous comparisons of these three modelling frameworks showed that all 
three produce identical results when the rates of cell proliferation and cell 
death arc sufficiently small relative to the rate of cell motility. Under these 
conditions the growth-to-conflucncc process takes place without developing 
any significant spatial correlations. Alternatively, under conditions where the 
proliferation and death rates are sufficiently large relative to the rate of cell 
motility, the growth-to-confluence process involves significant spatial correla- 
tions and the three models can make very different predictions owing to the 
extent to which each model includes a description of the effects of spatial cor- 
relations [1]. The key difficulty identified in our previous work was that each 
of the three models can always be calibrated to averaged experimental den- 
sity data to produce an estimate of the underlying rates of proliferation and 
death [1, 17]. This can be problematic since the calibration process always leads 
to a model prediction that matches the observed data, yet if the mean-field 
or moment-dynamics models are applied under inappropriate circumstances it 
is possible that the parameter estimates derived from them are meaningless. 
To address this issue we suggested that some measurement of the degree of 
spatial correlation ought to be considered to help make an informed decision 
about when different models ought to be applied [18,23]. 

Motivated by the importance of co-culture experiments, such as those 
shown in Figure l(d)-(f), the present work seeks to extend and generalize 
our previous study, which focused exclusively on mono-culture growth-to- 
confluence experiments, to now investigate appropriate modelling frameworks 
for analyzing co-culture growth-to-confluence experiments. We achieve this 
generalization by focusing on situations where we consider co-culture experi- 
ments with two cell types, which we refer to as cell type A and cell type B. 
This means that we always consider a total population composed of two sub- 
populations, subpopulation A and subpopulation B, and we anticipate that 
the general results and conclusions outlined here will also hold for more gen- 
eral cases involving three or more interacting subpopulations. In Section 2 we 
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outline three mathematical descriptions of a growth-to-confluence experiment 
with two cell types [11], and we show, in Section 3, that all three descriptions 
produce similar results when the rates of proliferation and death are suffi- 
ciently small whereas the three models produce very different results when 
the rates of proliferation and death are sufficiently large. These comparisons 
confirm that a simple model calibration procedure could lead to misleading 
results and in this light we suggest how additional measurements can be made 
to provide insight into how the most appropriate modelling framework could 
be chosen for a particular condition. We conclude with a brief discussion of 
our results in Section 4. 

2 Modelling methods 

In this work we will use three distinct mathematical models to describe the co- 
culture experiments: (i) a discrete model that explicitly incorporates individual 
cell behaviour; (ii) a traditional mean-field model that neglects spatial corre- 
lations; and (iii) a moment-dynamics model that approximately incorporates 
the effects of spatial correlations. Each of these models has been described in 
detail previously [11,15] and so we only provide a brief description of the key 
features of each of these models here. 

2.1 Discrete stochastic description 

We consider a population composed of two, possibly distinct, subpopulations 
on a two-dimensional square lattice, with lattice spacing A, and we invoke an 
exclusion mechanism whereby each lattice site can be occupied by, at most, 
a single agent [19]. Individual agents in each subpopulation undergo unbiased 
nearest neighbour motility, proliferation and death events. The i th species has 
a movement rate per unit time of P^, a proliferation rate per unit time of 
Pp and a death rate of P % d per unit time. To be consistent with a typical ex- 
perimental scenario, the lattice is initially uniformly populated, at random, 
meaning that site occupancies are initially uncorrelated and the density is, on 
average, spatially uniform. When an agent moves or proliferates, the target 
site is chosen at random from the relevant von Neumann neighbourhood, and 
the event is aborted if the target site is occupied. We invoke the simplest pos- 
sible proliferation mechanism whereby agents proliferate to form an identical 
daughter agent [19]. Simulations are propagated in time using a modified form 
of the Gillespie algorithm [6] as outlined in [11]. 

We report results from the stochastic model in two ways. First, we present 
visual snapshots showing the locations of agents in the population at different 
points in time. Second, we compute the average agent density across the lattice 
in the following way. If we are working on an L x x L y lattice, then we compute 




(1) 
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where Q\ (t) is the number of agents present from subpopulation i during 
the j th identically prepared realization of the same stochastic process, and 
we consider an ensemble of M realizations. In this way, (cj(t)) describes the 
average density of the i th subpopulation at time t. Since we focus on co-culture 
experiments with two different cell types, here i € [A, B]. Sample results from 
the stochastic model are shown in Figure 2 for a single-species model and in 
Figure 5 for a two-species model. Throughout this work we take L x = L y = 100 
and M = 100. 



2.2 Mean-field description 

The mean-field description of how the density of each population evolves over 
time can be written as 

- Pjcf, (2) 

where c™ f is the density of the i th subpopulation and i <G [A, B]. The super- 
script "mf explicitly refers to the fact that this model is based on making the 
traditional mean-field assumption whereby spatial correlations are completely 
neglected. The numerical solution of this system of two coupled ordinary dif- 
ferential equations is found using a fourth-order Runge-Kutta method with a 
constant time step [3]. 



dc 



ml' 



df 



P l c 



ml' 



,mf 



2.3 Moment-dynamics description 

The complete details of the derivation of the moment-dynamics description 
of this system was given previously by us [11] and so here we simply state 
the main results. The moment-dynamics model describing how the population 
evolves over time can be written 

- nc d , (3) 

where c™ f is the density of the i th subpopulation and i £ [A, B]. The super- 
script "md" explicitly refers to the fact that this model is based on an approx- 
imate moment-dynamics assumption whereby spatial correlations are approxi- 
mately incorporated using a moment closure assumption [11]. The i*y (Z\) term 
is a correlation function describing the correlation in occupancy of lattice sites 
at a distance of A between agents of subpopulation i and agents of subpop- 
ulation j. We note that if Fij(A) = 1 then the moment-dynamics description 
is equivalent to the mean-field description. For a two-species system with the 
total population being composed of subpopulation A and subpopulation B, 
this framework allow us to describe both the autocorrelation for each species, 



def 11 



= P'C 



md 
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Faa(A) and Fab(A), as well as the cross-correlation function, Fab(A), which 
by symmetry, is equivalent to Fba(A). 

To solve the moment-dynamics model we require equations governing the 
evolution of Fij(A), and these are generated by considering conservation equa- 
tions describing the evolution of finding pairs of agents separated by different 
lattice distances. In turn, these equations require description of the evolution 
of finding triplets of agents separated by different lattice distances and so on. 
This means that we have an infinitely large system of conservation equations 
which we must truncate to obtain an approximate solution. To this end, we 
close the system at second order using the Kirkwood superposition approxi- 
mation (KSA) [20]. 

On an infinite lattice, the partial differential equations describing evolution 
of the correlation functions are valid on the domain A < s < oo and are given 

by 



at 



1 _ V^ 2 „md 
1 l^k=l c k 



V 2 F«Ai(«) + F u {s) J2 4 nd V 2 F 4fe ( s ) 



fe=i 



X^X^Fuis) - 2P;X i {A)F ii (s), 



(4) 



and, for i ^ j, 
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where 

Y7 2 M = 

s ds V 9s 

The far-field boundary condition is 



Fy(s ->• oo) = 1, 



(7) 
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whereas the nearest-neighbour boundary conditions are given by 
dFu(A) P^ 



dt 



and, for i ^ j, 
dt 



Xi(A) 
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(9) 



with Xi(s) = 1 - Efc=i c k Fik(s). 

The numerical solution of this coupled system of nonlinear ordinary and 
partial differential equations is found by replacing the spatial derivative terms 
with a central difference approximation on a uniformly discretized domain with 
mesh size Sr. The resulting system of coupled nonlinear ordinary differential 
equations is integrated through time by discretising using the backwards Eulcr 
approximation, with time step St, and solving the resulting system of nonlinear 



s 
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algebraic equations using the Thomas algorithm [12] with Picard linearization 
and an absolute convergence tolerance of e [28]. 



3 Results 

We will present the results of our study in two sections: first, in Section 3.1, we 
present results for a single species problem which is relevant when considering 
a mono-culture growth-to-confluence experiment (Figure l(a)-(c)). Second, in 
Section 3.2, we present results for a two species problem which is relevant when 
considering a co-culture growth-to-confluence experiment (Figure l(d)-(f)). 



3.1 Single species results 

Results in Figure 2(a)-(c) compare the evolution of the average density profiles 
predicted by performing repeated stochastic simulations, and averaging the 
results, compared to the predictions of the mean-field and moment-dynamics 
models for a single species problem where the initial density of agents on 
the lattice is 10%. Results in Figure 2(a), corresponding to a high motility 
rate relative to the proliferation and death rates, show that both the mean- 
field and moment-dynamics models provide an excellent description of the 
averaged stochastic data. However, the results in Figure 2(b), corresponding 
to slightly lower motility rate, show that the mean-field model fails to ac- 
curately describe the evolution of the averaged stochastic data, whereas the 
moment-dynamics model produces results that are visually indistinguishable 
from the averaged stochastic data at this scale. These results indicate that 
the additional detail incorporated into the moment-dynamics model lead to 
an improved prediction under these circumstances. In contrast, the results in 
Figure 2(c), corresponding to zero motility, indicate that both the mean-field 
model and the more sophisticated moment-dynamics model can fail to pre- 
dict the evolution of the average density information. The reason why the 
mean-field and moment-dynamics models provide different results is because 
these two model frameworks make different assumptions about the underlying 
stochastic process. For example, the mean-field description completely neglects 
spatial correlations whereas the moment-dynamics model approximately ne- 
glects spatial correlations, and the accuracy of both descriptions decreases as 
the rate of proliferation and death increases relative to the rate of migration [1, 
16]. 

A key issue that arises when applying a mathematical model to interpret 
growth-to-conflucncc experiments, such as the images in Figure l(a)-(c), is 
that these images are often converted into a measure of average cell density as 
a function of time so that data is of the form presented in Figure 2(a)-(c) [2, 
24] . To estimate the proliferation and death rates of the cell population it then 
seems reasonable to calibrate these parameters using a relevant mathematical 
model so that the model predictions match the observed data. This kind of 
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Fig. 2 Results for the single species problem as the movement rate, P m , is varied. P m 
decreases from left to right: in the left-hand column P m = 500; in the centre column P m = 5; 
and in the right-hand column P m = 0. (a)— (c) Comparison of the averaged stochastic results 
(black) with predictions of the mean-field model (red) and the moment-dynamics model 
(green) over time, (d)-(f ) Snapshots from a single realisation of the stochastic model at t = 4. 
(g)— (i) Comparison of the averaged agent coordination number (red bars) with the predicted 
agent coordination number for a spatial uniform system (black asterisks). Parameters are 
as follows: P p = 1.0, P<j = 0.1 and initially 10% of lattice sites were occupied. 



calibration approach has been taken by many previous researchers in different 
contexts including the study of wound healing and malignant spreading [13, 
14,21]. One problem, highlighted by us previously [1,17], is that if we take 
average density information such as that in Figure 2(a)— (c), we can always 
calibrate either the mean-field or the moment-dynamics model to this data to 
produce estimates of P p and Pa- However, this commonly-invoked calibration 
process does not guarantee that the calibrated parameter estimates represent 
the actual birth and death rates since it is unclear, simply by inspecting density 
data alone, whether the mean- field and/or the moment-dynamics descriptions 
are valid. To provide insight into this question of model suitability we need to 
consider some additional information. For example, the snapshots of the dis- 
crete process, shown in Figure 2(d)-(f), illustrate how the agents are arranged 
on the domain at t = 4 for the parameter combinations associated with the 
density information in Figure 2(a)-(c), respectively. Visual inspection of these 
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images suggests that the agents are arranged relatively uniformly in Figure 
2(d), but that there is an increasing degree of spatial clustering and spatial 
correlation as the birth and death rates increase relative to the motility rate 
(Figure 2(e)-(f)). To help quantify these differences in spatial clustering and 
correlation we will now investigate two different quantities. 

3.1.1 Agent coordination number 

To provide a relatively straightforward measure of the degree of clustering and 
correlation in the distribution of agents at a particular time we will consider 
how the agent coordination number distribution, K, varies under different con- 
ditions. We define the agent coordination number, K, for a given site to be 
the total number of the eight closest sites in the Moore neighbourhood that 
are occupied, giving K £ [0,8]. For a randomly distributed, spatially uniform 
population at density C(t) £ [0,1], without any clustering or spatial corre- 
lations, the average agent coordination number at time t will be binomially 
distributed, 



with a mean of 8C(i) and a variance of 8C(t)(l — C(t)) [18]. This very straight- 
forward expression for the expected distribution of K could be compared 
with appropriately discretized images from an experiment, or from simula- 
tion data. The results in Figure 2(g)— (i) show, as a series of histograms, the 
agent coordination number distribution for the systems at t = 4 shown in Fig- 
ure 2(d)-(f). Superimposed on these histograms is the average observed agent 
coordination number. Comparing the results in Figure 2(g)— (i) we see that: 
(i) the observed agent coordination number compares well with the binomial 
distribution, equation (10), when motility is sufficiently high that the mean- 
field model accurately predicts evolution of the density (Figure 2(a)); (ii) we 
observe a relatively small discrepancy between the observed agent coordina- 
tion number and the binomial distribution when the motility is such that the 
moment-dynamics model captures the observed behaviour but the mean-field 
model fails (Figure 2(b)); and (iii) we observe a relatively large discrepancy 
between the observed agent coordination number and the binomial distribu- 
tion when motility is such that both the moment-dynamics model and the 
mean-field model fail to capture the observed behaviour (Figure 2(c)). 

The qualitative trends indicated in Figure 2 motivate us to consider quan- 
tifying these differences. First we define a measure of the difference between 
the uniform distribution of agent coordination number and the observed dis- 
tribution of agent coordination number 




(10) 



8 




(11) 



fc=0 



where ¥(K — k, t) is the averaged observed coordination number distribution 
at time t. To gain an understanding of how the difference in agent coordination 
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Fig. 3 A comparison of the errors between model and data for the single-species system. 
(a),(b) Evolution of Ej(t) for the mean- field model and the averaged stochastic data (j = mf, 
red), and for the moment-dynamics model and the averaged stochastic data (j = md, green) 
together with evolution of E^lt) (black). Parameters are as in the middle and left-hand 
columns of Figure 2, respectively, (c) The maximum value of Ej(t), Ej max , plotted as a 
function of the maximum value of Ek(t), ^kmaxi f° r j = m f (red) and j = md (green). The 
parameter ranges analysed axe Prn G [0, 500], P p G [0, 10], P d G [0, P p ] and the results from 
simulations using 140 different parameter combinations are displayed. The dashed horizontal 
lines represent 1% and 5% discrepancy in the density data. 



number is related to the difference in agent density, we also define 

E j (t) = \(c(t))-ci(t)\, (12) 

where j = mf when we are comparing the performance of the mean-field model, 
j = md when we are comparing the performance of the moment-dynamics 
model and (c(t)) is the averaged observed density at time t. 

One of the limitations of the results we presented in Figure 2(g)— (i) is that 
we compare coordination number data at one time point only. For completeness 
we present additional data, in Figure 3(a)-(b), where we consider single species 
growth-to-conflucncc experiments, with the same parameters used in Figure 
2(b)— (c), respectively, except now we provide data describing how £fc(t) and 
Ej{t) evolve with time for both the mean-field and moment-dynamics models. 
Results in Figure 3(a)-(b) illustrate some key features. Most notably, Ek(t) 
and Ej (t) vary dramatically with time, and reach some maximum value during 
the growth-to-confluence process. This means that we ought to be careful as 
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to the particular time we choose to measure Ek(t) and Ej(t), and we note that 
the maximum value of £7- (t) appears to occur earlier than the maximum value 
of Ej(t) in these cases (and all others we investigated). 

A summary of results in Figure 3(c) confirms that there is a clear trend 
between the maximum value of E).(t), -Ekmax, and the maximum value of 
Ej(t), -Ejmax, for both the mean- field and the moment-dynamics models. The 
deviation between the averaged discrete results and both the mean-field and 
moment-dynamics models increases as -Ekmax increases. Therefore, it is feasible 
to use a measure of -Ekmax to discriminate between the applicability of the 
mean-field, moment-dynamics and stochastic descriptions. The discrepancy 
for the mean-field model increases much more rapidly than the discrepancy 
for the moment-dynamics model. To illustrate how this kind of information 
might be used to distinguish between the suitability of the mean-field and 
moment-dynamics descriptions of this kind of process we have also included 
dashed horizontal lines in Figure 3(c) to show 1% and 5% thresholds in the 
discrepancy of the density data. For example, if we were content to use a model 
that provided no more than a 5% deviation between the averaged stochastic 
density data and the predictions of the model, then if the observed -Ekmax was 
less than approximately 0.25 a mean-field model would be suitable, whereas 
if the observed -Ekmax was between approximately 0.25 and 0.7, a moment- 
dynamics description would be suitable. If we observed -Ekmax > 0.7, then we 
could conclude that neither the mean-field or moment-dynamics descriptions 
were relevant and we ought to focus on using repeated stochastic simulation 
to interrogate our experimental data. 

3.1.2 Spatial correlation index 

As a second measure of the degree of clustering, we explore how the average 
correlations in lattice site occupancy for various lattice site distances vary with 
time. To obtain Faa(%> t) we count the number of pairs of agents separated by 
a distance x, and then divide this number by the square of the total number 
of agents on the lattice [10,23]: 

p / ^\ number of pairs of agents separated by distance x at time t 
' (number of agents at time t) 2 

(13) 

Figure 4(a) compares evolution of the correlation functions at distances A, 2 A 
and 3 A as a function of time against evolution of the error, Ej(t), between 
both the mean-field and moment-dynamics models and the averaged stochas- 
tic data. As with the agent coordination number, there is a delay between the 
maxima of the correlation function and the maximum error for both models. 
This indicates that we cannot rely on snapshots of the correlation functions 
to determine the suitability of a particular model. As for the agent coordi- 
nation number, we then explored how £j max depends on the maximum of 
FAA(x,t), -FAAmax, over a wide range of parameter values. Our results sug- 
gest that the values of i*AAmax and -Ej ma x are closely related for both models, 
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Fig. 4 Evolution of the correlation functions, (a) Evolution of the correlation functions 
at distances x = A, x = 2/A and x = 3A as a function of time (black lines), with the 
direction of increasing x indicated by the arrow. These results are superimposed on profiles 
showing Ej(t) for the mean-field model (red) and Ej(t) for the moment-dynamics model 
(green). Parameters are as in the middle column of Figure 2. (b)— (e) Plots of -Ej max as a 
function of the maximum of FAA{t), Fmax, with results from the mean- field model shown 
in red and from the moment dynamics model in green. The parameter ranges analysed are 
Pm 6 [0,500], P p e [0, 10], P d e [0,P P ] and the results from simulations using 90 different 
parameter combinations are displayed. 
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Fig. 5 Results for a two species problem, (a) The mean-field model (red) predicts both 
species evolve to occupy, on average, 50% of lattice sites, whilst the moment-dynamics 
model (green) correctly predicts the density of species A (solid line) to be greater than that 
of species B (dashed line). In each case the averaged discrete results (black) are obscured 
by the predictions of the moment-dynamics model. (b),(c) snapshots of the populations 
at t = 3 and t = 10, respectively, indicate that each of the species displays significant 
clustering. (d),(e) The averaged observed agent coordination number for each species at 
time t = 3. (f) The averaged auto- and cross-correlation functions at t = 3. The solid curves 
correspond to the solution of the moment-dynamics model whereas the dots correspond 
to averaged simulation data from the stochastic description. Parameters are as follows: 
= 20, PB = 5, P£ = 1, PB = 1, P£ = 0, Pf = 0 and each species initially occupies 
5% of lattice sites uniformly at random. 



over distances A, 2A and 3Z\, with the relationship becoming less well-defined 
as the distance increases (Figure 4(b)-(d)). The strongest trend corresponds 
to max{F J 4 J 4(2Z\) — Faa(A)} and Ej max (Figure 4(e)). Therefore, as an al- 
ternative to the agent coordination number data, described in Section 3.2.1, 
we suggest that estimates of max.{FAA(2A) — Faa(A)} could be used as a 
measure to distinguish between the suitability of the mean-field and moment- 
dynamics models in recapitulating the averaged discrete data for mono-culture 
growth-to-confluence experiments. 



3.2 Two species results 

We now focus on two species systems, which are relevant to co-culture growth- 
to-confluence experiments, like that shown in Figure l(d)-(f). Figure 5(a)-(c) 
shows results corresponding to simulations of a two species system, where the 
initial density of each species is 5%, so that the total population initially occu- 
pies 10% of the lattice, just like in Section 3.1 for the single species problem. 
Although both species A and B have the same proliferation rates, the evolution 
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of the two subpopulations is very different with subpopulation A, which has 
a higher movement rate compared to subpopulation B, occupying a greater 
percentage of lattice sites at later times. The evolution of the average density 
profiles in Figure 5(a) indicates that the mean-field model incorrectly predicts 
that the two subpopulations will eventually co-exist with each subpopulation 
occupying 50% of the lattice sites, whilst the moment-dynamics model cor- 
rectly predicts both the transient and steady state behaviour of the system 
in which subpopulation A dominates. The snapshots shown in Figure 5(b)- 
(c) of a single realisation at t = 3 and t — 10 indicates that the failure of 
the mean-field model may be due to the emergence of spatial correlations in 
lattice site occupancy, as was the case for the single-species model, since the 
build up of clustering in both subpopulations is visually distinct. Further in- 
vestigation into the failures of the mean-field and moment-dynamics models 
in various regions of parameter space indicate, as for the single-species model, 
that the validity of the model assumptions (the independency of lattice site 
occupancies and the validity of the KSA, respectively) [11]. 

As with the single species system, it is possible to calibrate either the mean- 
field or the moment-dynamics descriptions to a given set of averaged stochastic 
density data to provide an estimate of the model parameters governing the 
birth and death rates of both subpopulations in this context. However, as 
before, merely fitting different models may give rise to different parameter 
estimates, and it is unclear, a priori, which modelling framework is the most 
appropriate when we are dealing only with averaged density information, as is 
often the case in experimental cell biology [2,24]. To provide further guidance, 
we now investigate whether it is possible to use multi-species equivalents of the 
agent coordination number and spatial correlation index summary statistics 
to distinguish between different regions of parameter space in which the mean- 
field and moment-dynamics models are able to reliably replicate the averaged 
stochastic data. 



3.2.1 Agent coordination number 



xil-CiW-Cjit))*-*- 1 *. (14) 

This means that the agent coordination number of each individual subpopu- 
lation is binomially distributed, 



F(K t = ki,t) = Q)aw fc4 (i - c t {t)f~ k \ 



(15) 



as in the single species case. Results in Figure 5(d)-(e) show the observed 
averaged agent coordination numbers for the system shown in Figure 5(a)-(c) 
compared with the theoretical results for the case where the subpopulations 
are randomly distributed. Note that, had the mean-field model predicted the 
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Fig. 6 Averaged observed agent coordination number for two co-culture experiments, (a)— 
(c) In isolation, the agent coordination number of species A is close to that theoretically 
predicted for a uniform population, but this is not the case in co-culture. Parameters are as 



follows: P„ 



500, P£ 



5, Pi 



1 P h 

5 p 



1 P / 



: 0.1, 



Pi 



■ 0.1, t = 5 and both species 



initially occupy 5% of lattice sites uniformly at random, (d)— (f) Each individual species 
displays an agent coordination number that deviates significantly from that predicted for a 
uniform population, but this is not the case when they are considered as a single species. 



Parameters are as follows: P^ 



20, PJ; 



5, Pi 



1, Pi 



1 



pA 
^d 



- p ^' - 1 p -"-'-* d " ,J -' d 
and both species initially occupy 5% of lattice sites uniformly at random. 



0.1, Pf = 0.1, t 



observed average agent density data in Figure 5(a) correctly, then we would 
expect to see the actual distribution of agent coordination number match 
with the expected binomial result, similar to what we saw in Figure 2(g). In 
contrast, results in Figure 5(d)-(e) indicate that the observed distribution of 
agent coordination number deviates significantly from the binomial distribu- 
tion, which is consistent with the fact that the assumptions underlying the 
mean-field model are not met. 

At first glance, the results in Figure 5 indicate that the agent coordination 
number may be able to provide a guide as to the suitability of the mean- 
field and moment-dynamics models in predicting results from the stochastic 
model in the two-species co-culture model. As a potential means to simplify 
our calculations, we first ask whether understanding the dynamics of each 
subpopulation under equivalent mono-culture conditions is sufficient to deter- 
mine the ability of each modelling framework to replicate averaged stochastic 
data. Results in Figure 6 compare the observed and predicted distributions 
in agent coordination number using two different parameter combinations. In 
Figure 6(a)-(c), subpopulations A and B have the same parameter values as in 
the left-hand and centre columns of Figure 2, respectively. This means that, in 
isolation, we expect the evolution of population A to be well approximated by 
both the mean-field and moment-dynamics models, and that subpopulation B 
will be well-approximated by the moment-dynamics model and poorly approx- 
imated by the mean-field model. However, when these two subpopulations are 
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Fig. 7 A plot of -Ej max as a function of -Bjjmax over a wide range of parameter space reveals 
that the agent coordination number is not a suitable metric for distinguishing between 
the suitability of mean-field and moment-dynamics models in a co-culture experiment. The 
parameter ranges analysed axe Pfn £ [0,1000], P p e [0, 10], P d = 0 and the results from 
simulations using 55 different parameter combinations are displayed. 



grown together in co-culture, both subpopulations display agent coordination 
numbers that differ from that expected of a spatially uniform population even 
though the total agent coordination number indicates that, when considered 
as a single population, the system does not display significant clustering. This 
result implies that we cannot draw reliable conclusions regarding the suitabil- 
ity of a particular model for a co-culture experiment simply by considering the 
isolated behaviour of each of the individual species alone as in a mono-culture 
experiment. 

To provide further insight into modelling multi-species co-culture experi- 
ments we now explore whether it is necessary to be able to distinguish between 
the various subpopulations in order to determine the suitability of one mod- 
elling framework over another. Results in Figure 6(d)-(e) demonstrate that 
there are cases in which the distribution of agent coordination number for 
each subpopulation differ significantly from that expected from spatially uni- 
form populations. However, in Figure 6(f) we see that if we were unable to 
distinguish between the subpopulations, and instead we could only determine 
the total cell density, then calculating the agent coordination number for the 
total population would erroneously imply that the mean-field modelling frame- 
work could be appropriate. Together, these results suggest that, in order to 
use information about agent coordination number to distinguish between the 
suitability of the mean-field and moment-dynamics models, it is necessary to 
measure the agent coordination number of each subpopulation in a co-culture 
experiment. To deal with this complication, we would aim to produce a plot 
similar to that shown in Figure 3(c) where we showed that a strong correlation 
exists between differences in the predicted and observed densities and agent 
coordination numbers, £j ma x and -Ekmax, respectively, in the co-culture simu- 
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Fig. 8 A plot of -Ej max as a function of the maxima of the auto- and cross-correlation 
functions over a wide range of parameter space reveals that the agent coordination number 
is not a suitable metric for distinguishing between the suitability of mean-field and moment- 
dynamics models in a co-culture experiment. The parameter ranges analysed are Pm £ 
[0, 1000], P p £ [0, 10], Pd = 0 and the results from simulations using 45 different parameter 
combinations are displayed. 



lations. If this were possible, we could then reliably assume that measurements 
of -Ekmax for a particular co-culture system would be sufficient to guide our 
choice of model. Unfortunately, results shown in Figure 7 indicate that the re- 
sults are more subtle for the co-culture system: we see a much less pronounced 
relationship between £j ma x and -Ekmax m the co-culture system compared to 
the mono-culture system. Therefore, it does not seem reasonable to use the 
observed average agent coordination number as a means to distinguish the 
suitability of the mean-field and moment-dynamics models in the context of 
modelling a co-culture experiment. We now turn to study the spatial correla- 
tion index as an alternative summary statistic for co-culture experiments. 

3.2.2 Spatial correlation index 

Finally, we now explore whether the averaged observed correlation functions 
can improve our ability to make a distinction between the performance of the 
three modelling frameworks. To estimate Fij(x,t) we count the number of 
(i, j) agent pairs separated by a distance x, and normalise by the products of 
the total numbers of i and j agents on the lattice [10,23]: 

p r number of (?, j) agent pairs separated by a distance x at time t 

(number of i agents at time t) (number of j agents at time t) 

(16) 

In Figure 5(f) we demonstrate the ability of our moment-dynamics model 
to approximate both auto- and cross-correlations. As expected, the system 
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shows significant positive local spatial correlations in the distributions of each 
individual species, and negative local cross-correlations. In Figure 8 we plot 
-Ejmax as a function of the maxima of the auto- and cross-correlation functions, 
for a wide range of parameter values. As with the agent coordination number, 
the results are disappointing in that it does not seem sensible to based our 
model choice upon this measure of the system correlations in a co-culture 
experiment. 



4 Discussion and Conclusion 

In this work we explore the performance of three different mathematical mod- 
elling frameworks for quantitatively evaluating the results of in vitro mono- 
and co-culture growth-to-confluence assays. In particular, we focus on: (i) a 
stochastic description of a birth, death, movement process which is analyt- 
ically intractable but directly incorporates spatial correlation effects, (ii) a 
traditional mean-field description of a birth-death-movement process which 
is analytically tractable but completely neglects to incorporate any spatial 
correlation effects; and (iii) a more sophisticated moment-dynamics descrip- 
tion of a birth-death-movement process which is computationally tractable 
and approximately incorporates the effects of spatial correlations. Since it is 
relatively common to calibrate such models to experimental data in order to 
estimate model parameters [13,14,21], our aim is to provide a quantitative 
method which can be employed to distinguish between the validity of each 
type of model over a wide range of parameter space. 

For mono-culture assays in which we consider just one type of cell, we 
show that both the averaged agent coordination number and spatial correla- 
tion index provide a suitable means of distinguishing between the ability of 
each model to replicate the observed averaged density data over a wide range 
of parameter space. While our results for the mono-culture case suggest that 
we exercise caution against using just one time point during the experiment to 
measure the agent coordination number or spatial correlation index, we show 
that either of these spatial summary statistics can be used to reliably distin- 
guish between the application of a mean-field, moment-dynamics or stochastic 
representation of the growth-to-confluence process when we consider exam- 
ining the entire time course of the experiment and focus on the maximum 
deviation between the observed coordination number of spatial correlation in- 
dex for the experiment compared to the expected result for a system without 
any spatial correlations. 

In contrast, our results for co-culture assays suggest that great caution 
is warranted when calibrating mathematical models to replicate co-culture 
growth-to-confluence assays. We demonstrate that the spatial summary statis- 
tics for two species, applied in isolation, does not reliably indicate the suitabil- 
ity of a particular modelling framework when the two species are present in 
co-culture. Furthermore, we also show that extensions of the averaged agent 
coordination number data and spatial correlation index fail to provide a clear 
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quantitative measure of the suitability of one model over another for a wide 
range of parameter space. In this regard, our results suggest that extreme care 
ought be exercised when interpreting co-culture experiments using mathemat- 
ical models. Indeed, the results of this study suggest that further investigation 
of additional summary statistics, which can reliably distinguish between the 
application of stochastic, mean-field and moment-dynamics modelling frame- 
works, is warranted. 

In summary, our analysis indicates that it is possible to make a sensi- 
ble distinction between the suitability of mean-field, moment-dynamics and 
stochastic descriptions of a growth-to-confluence assay for a mono-culture 
growth-to-confluence experiment by analyzing either the distribution of agent 
coordination number or the spatial correlation index. Conversely, for a co- 
culture growth-to-confluence experiment more care is needed, and neither of 
these summary statistics can reliably distinguish between the suitability of the 
three candidate modelling frameworks. To this end, we suggest that the most 
pragmatic approach to make an initial assessment of the parameters governing 
a co-culture experiment is to use the discrete model as a screening tool to ob- 
tain parameter estimates. Depending on these initial estimates, it may then be 
possible to implement either a mean-field or a moment-dynamics description 
of the system for further analysis, if required. 
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